
% 
% HHVR_clogit_3_types_DM0.m
% Explain the subroutine here
% Note that it is computing only 2 types
% But is the limit of the 3 types sequence as DM --> infinite


function  [loglik  ] = HHVR_clogit_3_types_DM0(bbeta )  
   
  
global DATA choices   y_short X     individuals
 
s=size(DATA,1)/choices ;
num1=zeros(s,1);num2=zeros(s,1);num3=zeros(s,1);num4=zeros(s,1);num5=zeros(s,1);
prob_1=zeros(s,1);prob_2=zeros(s,1);prob_3=zeros(s,1);prob_4=zeros(s,1);prob_5=zeros(s,1);
den= zeros(s,1); den_x1= zeros(s,1);den_x2= zeros(s,1);den_x3=zeros(s,1);

score = zeros(s,size(bbeta,2));
size(score);


 for i=1:(size(DATA,1))/choices
  
     num1 (i) =  exp(X(i*choices-4,1)*bbeta(5) + X(i*choices-4,2)*bbeta(6)      ) ;
     num2 (i) =  exp(X(i*choices-3,1)*bbeta(5) + X(i*choices-3,2)*bbeta(6)   + 1*bbeta(1)  ) ;
     num3 (i) =  exp(X(i*choices-2,1)*bbeta(5) + X(i*choices-2,2)*bbeta(6)   + 1*bbeta(2)  ) ;
     num4 (i) =  exp(X(i*choices-1,1)*bbeta(5) + X(i*choices-1,2)*bbeta(6)   + 1*bbeta(3)  ) ;
     num5 (i) =  exp(X(i*choices  ,1)*bbeta(5) + X(i*choices  ,2)*bbeta(6)   + 1*bbeta(4)  ) ;
     
       den(i)   =  num1 (i)  +  num2 (i) + num3 (i) + num4 (i)  +  num5 (i)   ;
       prob_1(i)     =  num1 (i) /  den(i)  ;
       prob_2(i)     =  num2 (i) /  den(i)  ;
       prob_3(i)     =  num3 (i) /  den(i)  ;
       prob_4(i)     =  num4 (i) /  den(i)  ;
       prob_5(i)     =  num5 (i) /  den(i)  ;
       
       
     den_x1(i)   =  num1(i)*X(i*choices-4,1)  +  num2(i)*X(i*choices-3,1) + num3(i)*X(i*choices-2,1)+ num4(i)*X(i*choices-1,1)+ num5(i)*X(i*choices ,1)  ;
     den_x2(i)   =  num1(i)*X(i*choices-4,2)  +  num2(i)*X(i*choices-3,2) + num3(i)*X(i*choices-2,2)+ num4(i)*X(i*choices-1,2)+ num5(i)*X(i*choices ,2)  ; 
end
% prob = [prob_1  prob_2 prob_3 prob_4  prob_5  ] ;
% den_x  = [den_x1 den_x2  ];
%  
% prob2 = reshape(prob_2,1,1,size(prob_2,1));
% prob3 = reshape(prob_3,1,1,size(prob_3,1));
% prob4 = reshape(prob_4,1,1,size(prob_4,1));
% prob5 = reshape(prob_5,1,1,size(prob_5,1));
%  
% prob = [prob2 prob3 prob4  prob5  ]  ;
 
 

 
% What would look like in long code, and
% 2 x_ij bbetas
% a constant
% 8 choices
%  if ywill_(i)==1
%  score(i,1) = -  prob_2(i)*1;
%  score(i,2) = -  prob_3(i)*1;
%  score(i,3) = -  prob_4(i)*1;
%  score(i,4) = -  prob_5(i)*1;
%  score(i,5) = -  prob_6(i)*1;
%  score(i,6) = -  prob_7(i)*1;
%  score(i,7) = -  prob_8(i)*1;
%  score(i,8) =   (X(i*choices-6,1)- (den_x1(i) / den(i))) ;  
%  score(i,9) =   (X(i*choices-6,2)- (den_x2(i) / den(i))) ;
%  end

for i=1:(size(DATA,1))/choices
     
  
              if y_short(i)==1
             
              score(i,1)=     - ( num2 (i)/ den(i) )*1;
              score(i,2)=     - ( num3 (i)/ den(i) )*1;
              score(i,3)=     - ( num4 (i)/ den(i) )*1;
              score(i,4)=     - ( num5 (i)/ den(i) )*1;
              score(i,5) =   X(i*choices-4,1)- (den_x1(i) / den(i));
              score(i,6) =   X(i*choices-4,2)- (den_x2(i) / den(i));
               elseif y_short(i)==2
                  
              score(i,1)=     (1- (num2 (i)/ den(i))    *1);
              score(i,2)=     - ( num3 (i)/ den(i) )*1;
              score(i,3)=     - ( num4 (i)/ den(i) )*1;
              score(i,4)=     - ( num5 (i)/ den(i) )*1;
              score(i,5) =   X(i*choices-3,1)- (den_x1(i) / den(i));
              score(i,6) =   X(i*choices-3,2)- (den_x2(i) / den(i));
               elseif y_short(i)==3
                  
              score(i,1)=     - ( num2 (i)/ den(i) )*1;
              score(i,2)=     (1- ( num3 (i)/ den(i) )*1);
              score(i,3)=     - ( num4 (i)/ den(i) )*1;
              score(i,4)=     - ( num5 (i)/ den(i) )*1;
              score(i,5) =   X(i*choices-2,1)- (den_x1(i) / den(i));
              score(i,6) =   X(i*choices-2,2)- (den_x2(i) / den(i));
                   elseif y_short(i)==4
                  
              score(i,1)=     - ( num2 (i)/ den(i) )*1;
              score(i,2)=     - ( num3 (i)/ den(i) )*1;
              score(i,3)=     (1- ( num4 (i)/ den(i) )*1);
              score(i,4)=     - ( num5 (i)/ den(i) )*1;
              score(i,5) =   X(i*choices-1,1)- (den_x1(i) / den(i));
              score(i,6) =   X(i*choices-1,2)- (den_x2(i) / den(i));
                    elseif y_short(i)==5
                  
              score(i,1)=     - ( num2 (i)/ den(i) )*1;
              score(i,2)=     - ( num3 (i)/ den(i) )*1;
              score(i,3)=     - ( num4 (i)/ den(i) )*1;
              score(i,4)=     (1- ( num5 (i)/ den(i) )*1)   ;
              score(i,5) =   X(i*choices ,1)- (den_x1(i) / den(i));
              score(i,6) =   X(i*choices ,2)- (den_x2(i) / den(i));
                      
             else
             disp 'errorRR'
             pause
             end
           
 end
        
  
 
sscore = sum(score) ;
loglik = [sscore ];
 
 